﻿//#define ALGLIB
#define IMPLICIT_SHALLOW

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Text;

namespace alib.Vectors
{
	public static class vec_ext
	{
		///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		/// 
		///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		public static Double DotProduct(this double[] v, double[] v2)
		{
			Double sxy = 0;
			for (int i = 0; i < v.Length; i++)
				sxy += v[i] * v2[i];
			return sxy;
		}
		public static void Normalize(this double[] v)
		{
			Double l = VectorLength(v);
			if (l > 0.0)
				for (int i = 0; i < v.Length; i++)
					v[i] /= l;
		}
		public static Double VectorLength(this double[] v)
		{
			Double s = 0;
			foreach (Double n in v)
				s += n * n;
			return System.Math.Sqrt(s);
		}

		public static double[] MeanVector(this IEnumerable<double[]> seq)
		{
			var e = seq.GetEnumerator();
			double[] ret = null;
			int c = 0;
			while (e.MoveNext())
			{
				var cur = e.Current;
				if (ret == null)
				{
					ret = new double[cur.Length];
					cur.CopyTo(ret, 0);
				}
				else
				{
					for (int i = 0; i < ret.Length; i++)
						ret[i] += cur[i];

				}
				c++;
			}
			for (int i = 0; i < c; i++)
				ret[i] /= c;
			return ret;
		}
	}
}

namespace alib.Matrix
{
	using Math = System.Math;
	using String = System.String;

	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	/// 
	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	public static class _matrix_ext
	{
		public static T[] Row<T>(this T[,] v, int i)
		{
			T[] row = new T[v.GetLength(1)];
			for (int j = 0; j < row.Length; j++)
				row[j] = v[i, j];
			return row;
		}
		public static T[] Column<T>(this T[,] v, int j)
		{
			T[] col = new T[v.GetLength(0)];
			for (int i = 0; i < col.Length; i++)
				col[j] = v[i, j];
			return col;
		}

		public static T[,] Trim<T>(this T[,] m, int r0, int c0)
		{
			if (m.GetLength(0) <= r0 && m.GetLength(1) <= c0)
				return m;
			T[,] _new = new T[r0, c0];
			for (int i = 0; i < r0; i++)
				for (int j = 0; j < c0; j++)
					_new[i, j] = m[i, j];
			return _new;
		}

		public static T[,] TakeKRows<T>(this T[,] m, int k) { return Trim(m, k, m.GetLength(1)); }
		public static T[,] TakeKColumns<T>(this T[,] m, int k) { return Trim(m, m.GetLength(0), k); }
		public static T[,] TakeKRowsAndColumns<T>(this T[,] m, int k) { return Trim(m, k, k); }

		public static void NormalizeColumns(this double[,] m)
		{
			int c_rows = m.GetLength(0);
			int c_cols = m.GetLength(1);
			double[] lengths = new double[c_cols];
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					lengths[j] += Math.Pow(m[i, j], 2);
			for (int j = 0; j < c_cols; j++)
				lengths[j] = Math.Sqrt(lengths[j]);

			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					m[i, j] /= lengths[j];
		}

		public static T[,] Copy<T>(this T[,] m)
		{
			int c_rows = m.GetLength(0);
			int c_cols = m.GetLength(1);
			T[,] ret = new T[c_rows, c_cols];
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					ret[i, j] = m[i, j];
			return ret;
		}
	};

	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	/// 
	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	[DebuggerDisplay("{_DebugString()}")]
	public class Matrix
	{
		/// <summary> shallow (usurping) copy </summary>
		public Matrix(double[,] m)
		{
			this._values = m;
			this.c_rows = m.GetLength(0);
			this.c_cols = m.GetLength(1);
		}

		public Matrix(int c_rows, int c_cols)
			: this(new double[c_rows, c_cols])
		{
		}

		/// <summary> deep copy constructor </summary>
		public Matrix(Matrix m)
			: this(new double[m.c_rows, m.c_cols])
		{
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					_values[i, j] = m[i, j];
		}

		public Matrix(double[] v)
			: this(v.Length, v.Length)
		{
			for (int i = 0; i < v.Length; i++)
				_values[i, i] = v[i];
		}

		public Matrix(double[] v, int c_take)
			: this(c_take, c_take)
		{
			c_take = Math.Min(c_take, v.Length);
			for (int i = 0; i < c_take; i++)
				_values[i, i] = v[i];
		}

		protected readonly double[,] _values;
		protected readonly int c_rows;
		protected readonly int c_cols;

		public double this[int row, int column]
		{
			get { return _values[row, column]; }
			set { _values[row, column] = value; }
		}

		public int RowCount { get { return c_rows; } }
		public int ColumnCount { get { return c_cols; } }

		public double[] Row(int i)
		{
			double[] v = new double[c_cols];
			for (int j = 0; j < c_cols; j++)
				v[j] = _values[i, j];
			return v;
		}

		public double[] Column(int j)
		{
			double[] v = new double[c_rows];
			for (int i = 0; i < c_rows; i++)
				v[i] = _values[i, j];
			return v;
		}

		public IEnumerable<double[]> Rows
		{
			get
			{
				//int x = _values.Cast<Double>().Count();	// can access all in [,] this way
				for (int i = 0; i < c_rows; i++)
					yield return Row(i);
			}
		}

		public IEnumerable<double[]> Columns
		{
			get
			{
				for (int i = 0; i < c_cols; i++)
					yield return Column(i);
			}
		}

		public void CleaveVertical(int left_cols, out Matrix left, out Matrix right)
		{
			if (left_cols == c_cols)
			{
				left = new Matrix(this);
				right = new Matrix(c_rows, 0);
			}
			else if (left_cols == 0)
			{
				left = new Matrix(c_rows, 0);
				right = new Matrix(this);
			}
			else if (left_cols > c_cols)
				throw new Exception();
			else
			{
				left = new double[c_rows, left_cols];
				right = new double[c_rows, c_cols - left_cols];
				for (int i = 0; i < c_rows; i++)
					for (int j = 0; j < c_cols; j++)
						if (j < left_cols)
							left[i, j] = _values[i, j];
						else
							right[i, j - left_cols] = _values[i, j];
			}
		}

		public void CleaveHorizontal(int top_rows, out Matrix top, out Matrix bottom)
		{
			if (top_rows == c_rows)
			{
				top = new Matrix(this);
				bottom = new Matrix(0, c_cols);
			}
			else if (top_rows == 0)
			{
				top = new Matrix(0, c_cols);
				bottom = new Matrix(this);
			}
			else if (top_rows > c_rows)
				throw new Exception();
			else
			{
				top = new double[top_rows, c_cols];
				bottom = new double[c_rows - top_rows, c_cols];
				for (int i = 0; i < c_rows; i++)
					for (int j = 0; j < c_cols; j++)
						if (i < top_rows)
							top[i, j] = _values[i, j];
						else
							bottom[i - top_rows, j] = _values[i, j];
			}
		}

		public void SetColumnValues(int j, double[] v)
		{
			if (v.Length != c_rows)
				throw new Exception();
			for (int i = 0; i < c_rows; i++)
				_values[i, j] = v[i];
		}

		public static Matrix Identity(int size)
		{
			Matrix m = new Matrix(size, size);
			for (int i = 0; i < size; i++)
				for (int j = 0; j < size; j++)
					m[i, j] = (i == j) ? 1.0 : 0.0;
			return m;
		}

		public Matrix Transpose()
		{
			Matrix m = new Matrix(c_cols, c_rows);
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					m[j, i] = this[i, j];
			return m;
		}

		public static Matrix Add(Matrix m_l, Matrix m_r)
		{
			Debug.Assert(m_l.c_cols == m_r.c_cols);
			Debug.Assert(m_l.c_rows == m_r.c_rows);

			Matrix m = new Matrix(m_l.c_rows, m_r.c_cols);
			for (int i = 0; i < m_l.c_rows; i++)
				for (int j = 0; j < m_l.c_cols; j++)
					m[i, j] = m_l[i, j] + m_r[i, j];
			return m;
		}

		public static Matrix Subtract(Matrix m_l, Matrix m_r)
		{
			Debug.Assert(m_l.c_cols == m_r.c_cols);
			Debug.Assert(m_l.c_rows == m_r.c_rows);
			Matrix m = new Matrix(m_l.c_rows, m_r.c_cols);
			for (int i = 0; i < m_l.c_rows; i++)
				for (int j = 0; j < m_l.c_cols; j++)
					m[i, j] = m_l[i, j] - m_r[i, j];
			return m;
		}

		public static Matrix Multiply(Matrix m_l, Matrix m_r)
		{
			if (m_l.c_cols != m_r.c_rows)
				throw new Exception();
			Matrix m = new Matrix(m_l.c_rows, m_r.c_cols);
			for (int i = 0; i < m.c_cols; i++)
			{
				for (int j = 0; j < m_l.c_rows; j++)
				{
					double value = 0.0;
					for (int k = 0; k < m_r.c_rows; k++)
						value += m_l[j, k] * m_r[k, i];
					m[j, i] = value;
				}
			}
			return m;
		}

		public static Matrix Multiply(double left, Matrix m_r)
		{
			Matrix m = new Matrix(m_r.c_rows, m_r.c_cols);
			for (int i = 0; i < m.c_rows; i++)
				for (int j = 0; j < m_r.c_cols; j++)
					m[i, j] = left * m_r[i, j];
			return m;
		}

		public static Matrix Multiply(Matrix m_l, double right)
		{
			Matrix m = new Matrix(m_l.c_rows, m_l.c_cols);
			for (int i = 0; i < m_l.c_rows; i++)
				for (int j = 0; j < m_l.c_cols; j++)
					m[i, j] = m_l[i, j] * right;
			return m;
		}

		public static Matrix Divide(Matrix m_l, double right)
		{
			Matrix m = new Matrix(m_l.c_rows, m_l.c_cols);
			for (int i = 0; i < m_l.c_rows; i++)
				for (int j = 0; j < m_l.c_cols; j++)
					m[i, j] = m_l[i, j] / right;
			return m;
		}

		public static double[] DiagonalVector(Matrix m)
		{
			Debug.Assert(m.c_rows == m.c_cols);

			double[] v = new double[m.c_rows];
			for (int i = 0; i < m.c_rows; i++)
				v[i] = m[i, i];
			return v;
		}

		public Matrix Trim(int r0, int c0)
		{
			if (c_rows <= r0 && c_cols <= c0)
				return new Matrix(this);
			double[,] _new = new double[r0, c0];
			for (int i = 0; i < r0; i++)
				for (int j = 0; j < c0; j++)
					_new[i, j] = _values[i, j];
			return new Matrix(_new);
		}

		public Matrix TakeKRows(int k) { return Trim(k, c_cols); }
		public Matrix TakeKColumns(int k) { return Trim(c_rows, k); }
		public Matrix TakeKRowsAndColumns(int k) { return Trim(k, k); }

		public void NormalizeColumns()
		{
			double[] lengths = new double[c_cols];
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					lengths[j] += Math.Pow(_values[i, j], 2);
			for (int j = 0; j < c_cols; j++)
				lengths[j] = Math.Sqrt(lengths[j]);

			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					_values[i, j] /= lengths[j];
		}
		public void NormalizeRows()
		{
			double[] lengths = new double[c_rows];
			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					lengths[i] += Math.Pow(_values[i, j], 2);
			for (int i = 0; i < c_rows; i++)
				lengths[i] = Math.Sqrt(lengths[i]);

			for (int i = 0; i < c_rows; i++)
				for (int j = 0; j < c_cols; j++)
					_values[i, j] /= lengths[i];
		}

		public static Matrix operator +(Matrix m_l, Matrix m_r) { return Matrix.Add(m_l, m_r); }
		public static Matrix operator -(Matrix m_l, Matrix m_r) { return Matrix.Subtract(m_l, m_r); }
		public static Matrix operator *(Matrix m_l, Matrix m_r) { return Matrix.Multiply(m_l, m_r); }
		public static Matrix operator *(double left, Matrix m_r) { return Matrix.Multiply(left, m_r); }
		public static Matrix operator *(Matrix m_l, double right) { return Matrix.Multiply(m_l, right); }
		public static Matrix operator /(Matrix m_l, double right) { return Matrix.Divide(m_l, right); }

#if IMPLICIT_SHALLOW
		public static implicit operator Matrix(double[,] m) { return new Matrix(m); }
		public static implicit operator double[,](Matrix m) { return m._values; }
#else
		public static implicit operator Matrix(double[,] m) { throw new NotImplementedException(); }
		public static implicit operator double[,](Matrix m)
		{
			double[,] result = new double[m.c_rows, m.c_cols];
			for (int i = 0; i < m.c_rows; i++)
				for (int j = 0; j < m.c_cols; j++)
					result[i, j] = m[i, j];
			return result;
		}
#endif

#if ALGLIB
		/// <summary>
		/// Make alglib's SVD nicer by using the above type
		/// </summary>
		public static bool SingularValueDecomposition(Matrix A, out double[] w, out Matrix V)
		{
			w = null;
			V = null;
			double[,] u = null;
			double[,] v = null;

			if (!alglib.svd.rmatrixsvd(A, A.c_rows, A.c_cols, 0, 1, 2, ref w, ref u, ref v))
				return false;
			V = new Matrix(v);
			return true;
		}
#endif

		public String TextDisplay()
		{
			StringBuilder sb = new StringBuilder();
			for (int j = 0; j < c_cols; j++)
				sb.AppendFormat("{0,6}", j);
			sb.AppendLine();
			for (int j = 0; j < c_cols; j++)
				sb.AppendFormat(" -----");
			sb.AppendLine();
			for (int i = 0; i < c_rows; i++)
			{
				for (int j = 0; j < c_cols; j++)
				{
					//String s = String.Format("{0,6:0.00}", _values[i, j]);
					//if (s.EndsWith("0.00"))
					//    s = "      ";
					//int x = (int)(_values[i, j] * 100.0);
					//String s = String.Format("{0,6}", x == 0 ? "" : x.ToString());
					//sb.Append(s);
					var d = _values[i, j];
					String s;
					if (alib.Math._math_ext.IsZero(d))
					{
						s = "";
					}
					else
					{
						s = String.Format("{0:N}", d);
					}
					sb.Append(alib.String._string_ext.TrimEndOrPadLeft(s, 6));
				}
				sb.AppendLine();
			}
			return sb.ToString();
		}

		private String _DebugString() { return String.Format("{0} rows by {1} columns", c_rows, c_cols); }

		public override String ToString()
		{
			return TextDisplay();
		}
	};
}

